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COMPUTATION  OF  NONLINEAR  BACKSCATTERING  USING  A  HIGH-ORDER 

NUMERICAL  METHOD* 

G.  FIBICHU,  B.  ILANt§?  AND  S.  TSYNKOV^ 

Abstract,  The  nonlinear  Schrodinger  equation  (NLS)  is  the  standard  model  for  propagation  of  intense 
laser  beams  in  Kerr  media.  The  NLS  is  derived  from  the  nonlinear  Helmholtz  equation  (NLH)  by  employing 
the  paraxial  approximation  and  neglecting  the  backscattered  waves.  In  this  study  we  use  a  fourth-order 
finite-difference  method  supplemented  by  special  two-way  artificial  boundary  conditions  (ABCs)  to  solve  the 
NLH  as  a  boundary  value  problem.  Our  numerical  methodology  allows  for  a  direct  comparison  of  the  NLH 
and  NLS  models  and  for  an  accurate  quantitative  assessment  of  the  backscattered  signal. 

Key  words,  Kerr  medium,  wave  propagation,  self  focusing,  fourth-order  method,  two-way  ABCs 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  The  propagation  of  intense  laser  beams  (time-harmonic  electromagnetic  waves)  in  a 
bulk  Kerr  medium  is  usually  modeled  by  the  critical  nonlinear  Schrodinger  equation  (NLS)  for  the  electric 
field  amplitude.  Since  light  rays  bend  toward  areas  with  higher  index  of  refraction,  the  nonlinear  dependence 
of  the  index  of  refraction  on  beam  intensity  has  a  self-focusing  effect,  whereby  a  sufficiently  intense  laser 
beam  becomes  narrower  as  it  propagates.  In  particular,  the  NLS  model  predicts  that  when  the  input  beam 
power  ( L 2  norm)  exceeds  a  given  critical  threshold  NC7  then  the  beam  can  collapse  to  a  point  at  a  finite 
propagation  distance.  For  more  information  of  self-focusing,  see  e.g.,  [6,9]. 

As  the  beam  propagates  it  induces  changes  in  the  optical  properties  of  the  medium.  As  a  result,  part  of 
the  incoming  wave  is  reflected  back,  a  phenomenon  referred  to  as  backs cattering.  Very  little  is  actually  known 
on  backscattering  in  nonlinear  self-focusing,  except  for  the  general  belief  that  it  is  “small.”  Since,  however, 
small-magnitude  mechanisms  can  have  a  large  effect  in  self-focusing  [6],  there  is  a  need  to  accurately  quantify 
the  magnitude  of  backscattering  and  study  how  this  phenomenon  may  affect  the  beam  propagation.  Another 
application  which  could  greatly  benefit  from  better  understanding  of  backscattering  is  remote  sensing  of  the 
atmosphere  [12],  where  the  measured  signal  is  exactly  the  backscattered  wave. 

The  backscattered  wave  is  neglected  in  the  NLS  model  which  only  describes  the  forward-propagating 
wave.  Calculation  of  backscattering  requires,  therefore,  going  back  to  the  nonlinear  Helmholtz  equation 
(NLH),  from  which  the  NLS  is  derived.  The  NLS  is  an  evolution  equation  with  the  spatial  coordinate  in  the 
direction  of  propagation  playing  the  role  of  “time.”  Therefore,  the  correct  mathematical  formulation  for  the 
NLS  is  the  Cauchy  (initial  value)  problem  and  as  such,  solving  it  numerically  is  a  relatively  straightforward 
computational  procedure.  In  contradistinction  to  that,  the  NLH  is  elliptic  in  its  nature,  and  a  special 
multidimensional  boundary  value  problem  needs  to  be  formulated  and  solved  for  this  equation,  which  is  a 
much  harder  task  from  the  standpoint  of  computing.  The  first  numerical  simulations  based  on  solving  a 
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boundary  value  problem  for  the  NLH  were  recently  performed  in  [7]  using  an  advanced  fourth-order  method. 
In  that  study,  the  design  fourth-order  convergence  rate  of  the  method  was  corroborated  experimentally  on  a 
model  linear  problem.  Subsequently,  a  series  of  the  grid  convergence  tests  were  conducted  in  the  nonlinear 
regime.  In  the  current  paper  we  go  beyond  grid- convergence  arguments  and  show  that  the  asymptotic  limit  of 
the  NLH  solutions  obtained  in  the  simulations  is  the  corresponding  NLS  solution.  This  comparison  provides 
strong  support  that  the  calculated  NLH  solution  is  indeed  the  physical  one.  In  Section  2.3  we  obtain  an 
asymptotic  estimate  of  the  magnitude  of  backscattering  and  subsequently  show  in  Section  4  that  it  agrees 
with  the  calculated  values. 

When  the  initial  datum  is  sufficiently  large  the  NLS  solution  develops  singularities  at  a  finite  propagation 
distance  (see  Section  2.5).  Since,  however,  physical  quantities  do  not  become  infinite,  a  natural  question 
is  whether  the  corresponding  solution  of  the  NLH  exists  globally.  This  fundamental  question  has  been 
open  for  many  years.  There  have  been  indications,  though,  that  solutions  to  the  NLH  exist  even  when 
the  corresponding  NLS  solutions  become  singular,  based  on  both  numerical  solution  of  “modified”  NLH 
equations  [1,2,4]  and  on  asymptotic  analysis  [5],  but  these  studies  did  not  take  into  account  backscattering 
effects.  Therefore,  our  long-term  goal  is  to  solve  the  NLH  for  those  incoming  signals  that  lead  to  blowup 
in  the  NLS  model.  In  the  current  study,  however,  we  concentrate  on  the  more  attainable  goal  of  better 
understanding  (in  terms  of  both  analysis  and  numerical  simulations)  the  regime  when  the  corresponding 
solution  of  the  NLS  does  not  blow  up.  Our  hope  is  that  this  understanding  will  eventually  allow  to  solve  the 
NLH  for  “any”  incoming  signal. 

2.  Mathematical  Models. 

2.1.  The  Nonlinear  Helmholtz  Equation.  A  typical  experimental  setup  (both  physical  and 
numerical)  for  the  propagation  of  waves  in  Kerr  media  is  shown  in  Figure  2.1.  An  incoming  laser  beam 
with  known  characteristics  impinges  normally  on  the  planar  interface  z  —  0  between  the  linear  and  the 
nonlinear  media.  The  electric  field  E  —  E(x±, . . .  ,£z>_i,  z)  in  WP  is  governed  by  the  nonlinear  Helmholtz 
equation 

(2.1)  (dzz  +  A±)E  +  k2E  =  0,  k2  =  k%(l  +  e|i?|2<r),  (an, . . . ,xD-i)  e  R0-1,  z>0, 

where  k o  is  the  wavenumber,  e  =  4tocri27  n 2  is  the  Kerr  coefficient,  and  Ax  =  dXlXl  +  ...  +  dXD  lXD  x  is 
the  transverse  Laplacian  (the  diffraction  term),  see,  e.g.,  [3,8].  For  simplicity  we  consider  from  now  on  the 
cylindrically-syrrirnetric  case  where  E  =  E(r7  z)  and  r  =  yjx\  +  . . .  +  xjj^ . 

The  nonlinear  medium  occupies  the  semi-space  z  >  0  (see  Figure  2.1).  Consequently,  the  NLH  (2.1) 
has  to  be  supplemented  by  boundary  conditions  at  z  =  0  and  z  — >  +00.  We  require  that  as  z  — >  +00,  E 
has  no  left-traveling  components  and  that  the  propagation  is  diffraction-dominated  with  the  field  amplitude 
decaying  to  zero,  i.e.,  lirn  max  \E(r?z)\  =  0,  which  also  means  lirri  k2  =  k%  .  In  other  words,  at  large  z7 s 

z—too  0<r<oo  z-y+oo 

the  solution  should  be  a  linear  superposition  of  right-traveling  waves.  Since  the  actual  numerical  simulation 
is  carried  out  on  a  truncated  domain  0  <  z  <  zm ax  (Figure  2.1),  the  desired  behavior  of  the  solution  as 
z  — y  +00  has  to  be  captured  by  a  far-field  artificial  boundary  condition  (ABC)  at  the  artificial  boundary 
z  =  zmax.  This  boundary  condition  should  guarantee  a  reflectionless  propagation  of  all  the  waves  traveling 
toward  z  =  +00.  Often,  boundary  conditions  designed  to  ensure  the  transparency  of  the  outer  boundary  to 
the  outgoing  waves  are  called  radiation  boundary  conditions  [10]. 

The  situation  is  more  complicated  at  the  interface  z  —  0,  where  the  total  field  E(r,  0)  is  composed 
of  a  given  incoming  (right-traveling)  component  Elnc(r70)  and  an  unknown  backscattered  (left-traveling) 
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Fig.  2.1.  Schematic  of  propagation  of  waves  in  Kerr  media. 


component  l?scat(r?0)?  i.e., 

E(r,  0)  =  Einc(r,  0)  +  .Escat(r,  0)  . 

As  such,  the  boundary  condition  at  z  =  0  has  to  guarantee  the  reflectionless  propagation  of  any  left- 
traveling  wave  through  the  interface  and  at  the  same  time  be  able  to  correctly  prescribe  the  incoming  signal. 
Implementation  of  such  two-way  ABC  was  done  in  [7]. 

Finally,  we  assume  symmetry  at  r  —  0  and  vanishing  of  the  electric  field  as  r  — >  +oo.  In  practice,  we 
truncate  the  domain  at  some  large  but  finite  rmax  and  require  that  E(rmaxyz)  =  0.  The  justification  for  the 
use  of  this  approach  to  treat  the  “lateral”  boundaries  can  be  found  in  [10]  and  the  bibliography  there. 

Let  us  also  note  that  in  this  study  we  do  not  take  into  account  the  effect  of  discontinuity  in  the  index 
of  refraction  across  the  interface  z  =  0.  Thus,  we  assume  that  Emc(r?  0)  is  the  incoming  wave  after  it  has 
already  passed  through  the  interface  (i.e.,  at  z  —  0+).  We  also  assume  that  left-traveling  (i.e.,  backscattered) 
waves  are  not  reflected  by  the  interface  z  =  0  back  into  the  domain  z  >  0.  One  can  expect  the  latter  effect 
to  be  small;  it  will  be  investigated  in  a  future  study. 

2.2.  Paraxial  Approximation  and  the  Nonlinear  Schrodinger  Equation.  Let  be  the  initial 
width  of  the  impinging  beam.  We  first  introduce  the  dimensionless  quantities  f,  z7  and  ift  as 

(2.2)  f  =  —  ,  E  =  eik°*(er*kt)-1/2*i>(r,z)  , 

ro  ILof 

where  L of  =  ^oro  is  the  diffraction  length .  Then,  by  dropping  the  tildes  we  obtain 

(2.3)  iipz  +  Aj.V’  +  =  -4 f2tpzz  , 

where  /  =  l/r^k q  1  is  the  nonparaxiality  parameter . 

The  standard  derivation  of  the  NLS  is  based  on  the  assumption  that  the  envelope  ^  is  slowly  varying. 
In  that  case,  one  can  neglect  the  term  on  the  right-hand  side  of  (2.3)  [i.e.,  apply  the  paraxial  approximation] 
and  obtain  the  nonlinear  Schrodinger  equation 

(2.4)  i'tftz  +  A +  \i>\2<Ti>  =  0  * 

The  NLS  (2.4)  is  supplemented  by  the  initial  condition  at  z  =  0 

i>{r,  0)  =  (er^)1/2<r£’inc(r,0)  . 
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Subsequently,  it  needs  to  be  integrated  by  a  “time” -marching  algorithm,  where  the  direction  of  propagation  z 
plays  the  role  of  time.  We  reemphasize  that  backs cattering  effects  are  not  taken  into  account  by  the  NLS  (2.4)* 
Indeed,  once  (2.4)  is  solved,  the  overall  solution,  according  to  (2.2),  is  the  slowly  varying  amplitude  times 
the  forward  propagating  oscillatory  component  elk°z . 

2.3.  Preliminary  Analysis  of  Backscattering,  To  the  best  of  our  knowledge,  no  accurate 
quantitative  analysis  of  backscattering  in  nonlinear  self-focusing  has  ever  been  performed,  although  there  is  a 
general  belief  that  the  magnitude  of  the  backscattered  signal  is  small.  In  this  section  we  present  a  preliminary 
asymptotic  study  of  backscattering.  To  do  that,  we  consider  a  more  general  ansatz  for  E  than  (2.2)  which 
is  composed  of  both  forward-propagating  and  backward- propagating  waves,  i.e., 


(2.5) 


E  =  (er2k 


A(r,  z)eikoZ  +  B(r,  z)e 


—ikoz 


where  A  and  B  are  slowly- varying  envelopes.  Substitution  in  the  NLH  (2.1)  yields 


Jkoz 


A„  +  2  ikoAz  +  Aj_A  +  \A  +  e~2ik°z  B\2*  A 


+  e 


—ikoz 


Bzz  +  2ik0Bz  +  A  ±B  +  \A  +  e~2ik°z  Bf  B 


0  . 


Changing  to  the  nondimensional  variables  (2.2)  gives  (after  dropping  the  tildes) 


JWf  > 


fAzz  +  iAz  +  Aj_A  +  \A  +  e~iW  f2)z. B\2<r  A  +  fBzz  +  iBz  +  A  ±B  +\A  +  e_i(4/  f2)z  Bf  B 


0  . 


Let  us  average  the  last  equation  over  one  fast  oscillation.  For  example,  using  Taylor  expansion,  we  obtain 


7T/2 


H-tt/74 


-nf2/ 4 


?+tt/2/4 


-tt/2/4 


[B(z)  +  (C  -  z)Bz(z)  +  . . .]  d(  =  B(z)(  1  +  0(/2))  . 


Similarly, 


7 r/2 


■+*fl  4 


-7T/V4 


A(C  )dC  =  ^ 


-+wf2/i 


eWfa)*[A(z)  +  (C  -  z)Az(z)  +  . . .} d( 


=  ^AzeWf2>  +  0(f)  . 


Consequently,  we  obtain  the  following  equation  for  the  backscattered  wave 


(2.6) 

where 


f2Bzz  +  iBz  +  A ±B  +  \A  +  =  f2F  ,  B(z) 


=  0 


F  —  — 
% 


f2Azz  +  iAz  +  Aj_A  +  \A  +  e-i{VC)*BfA 


Let  us  now  employ  the  common  assumption  that  backscattering  is  small,  i.e.,  B  <C  A,  Since  /  1,  equation 

(2.6)  for  B  can  be  approximated  with  the  linear  Schrodinger  equation 


iBz  +  A±B+\A\2*B  =  f2F,  B(z) 


=  0, 


where 


F  —  —  e*(4/72)* 
4  i 


iAz  +  A±A  +  |A|JcrA 
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Since  the  solution  of 


iBz  +  A±B+\A\2°B  =  0  ,  B(z) 
is  B  =  0,  the  above  analysis  suggests  that 
(2.7)  B/A  =  0(f2)  . 


=  0 


If  the  quadratic  scaling  law  (2.7)  can  be  confirmed  independently,  then  it  will  provide  a  convincing 
justification  for  the  assumption  that  backscattering  is  indeed  small.  Hence,  in  our  simulations  we  expect  to 
see  that 


(2.8) 


(ajgy^\E\  - \A\ 

\A\ 


0(f). 


The  numerical  results  of  Section  4  do  corroborate  these  expectations. 
The  above  analysis  also  shows  that 


\E\  =  (er2k2)  1/2'l^l 


—2ikoz 


Therefore,  the  amplitude  of  the  NLH  solution  may  have  0(f 2)  ripples  with  the  wavelength  rr/ko  due  to 
backscattering  on  top  of  the  slowly  varying  amplitude  of  the  forward-propagating  wave.  NLH  simulations 
suggesting  that  this  indeed  may  be  the  case  have  been  reported  in  [7]  (see  also  Figure  4.2  in  Section  4).  It 
is  not  clear,  however,  to  what  extent  the  ripples  observed  in  [7]  are  a  numerical  artifact  due  to  placing  the 
far-field  artificial  boundary  too  close  to  the  nonlinear  self-focusing  zone.  Numerical  study  conducted  in  [7] 
did,  in  fact,  involve  the  analysis  of  how  the  location  of  the  far-field  artificial  boundary  affects  the  solution; 
placing  this  boundary  further  and  further  away  caused  the  reduction  in  the  ripples’  magnitude  but  never 
allowed  to  eliminate  them  completely.  This  may  still  imply,  though,  that  the  boundary  was  not  “sufficiently 
far”  away.  Therefore,  no  definite  conclusion  as  to  the  presence  of  the  0(f2)  ripples  in  the  NLH  solutions 
can  be  made  at  this  time  and  this  question  requires  a  subsequent  thorough  study. 

2.4.  Nonpar axiality  and  backscattering.  The  traditional  way  of  introducing  the  paraxial 
approximation  is  reported  in  Section  2.2,  where  the  right-hand  side  of  equation  (2.3)  is  omitted  and  the 
NLS  is  derived.  The  more  careful  analysis  of  Section  2.3  shows,  however,  that  two  approximations  are,  in 
fact,  being  made  when  the  NLH  is  approximated  with  the  NLS:  Neglecting  Azz  (the  paraxial  approximation 
in  the  narrow  sense,  i.e.,  as  it  applies  to  the  forward-propagating  waves)  and  neglecting  B  (backscattering). 
We  recall  that  previous  studies  [4,5]  suggested  that  nonpar  axiality  of  the  right-traveling  waves  (i.e.,  Azz  in 
the  sense  of  Section  2.3)  arrests  the  collapse  of  the  NLS  solutions,  but  these  studies  did  not  take  into  account 
backscattering  effects.  Having  said  that,  we  still  note  that  the  separation  into  nonparaxial  and  backscattering 
effects,  which  is  based  on  the  arisatz  (2.5),  is  somewhat  artificial,  since  the  problem  in  nonlinear.  Therefore, 
when  we  compare  the  NLH  and  NLS  solutions,  it  is  not  precisely  clear  which  part  of  the  difference  conies  from 
nonparaxial  effects  for  the  right-traveling  waves,  and  which  one  from  backscattering.  A  notable  exception  is7 
however 7  at  z  —  07  where  the  difference  between  NLH  and  NLS  solutions  is  solely  due  to  backscattering . 

2.5.  Critical  NLS.  It  is  well  known  that  solutions  of  the  NLS  (2.4)  can  become  singular  when  either 
<j(D  —  1)  >  2,  the  supercritical  NLS,  or  when  cr(D  —  1)  =  2,  the  critical  NLS  (D  is  the  space  dimension). 
However,  whereas  in  supercritical  collapse  nonlinearity  dominates  over  diffraction  near  the  singularity,  in 
the  critical  collapse  nonlinearity  and  diffraction  are  almost  balanced  near  the  singularity.  Consequently,  the 
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singularity  formation  is  highly  sensitive  to  small  perturbations  in  the  critical  case,  but  much  less  so  is  the 
supercritical  case. 

The  physical  case  that  corresponds  to  the  propagation  of  laser  beams  in  bulk  Kerr  media  is  the  critical 
one,  as  D  =  3  and  a  =  1.  However,  in  order  to  reduce  the  complexity  of  the  computations,  below  we  consider 
the  critical  case  D  =  2  and  a  =  2.  Thus,  the  NLH  for  E  —  E(r7z)  and  the  NLS  for  ip  =  ip(r7  z)7  which  are 
solved  numerically  in  this  study  are 

(2.9)  EZz  +  Err  +  kS(l  +  e\E\i)E  =  0, 
and 

(2.10)  ifz  +  lprr  +  |^>| V  =  0  j 

respectively. 

3.  Numerical  Methods.  The  NLH  (2.9)  is  solved  using  fourth-order  finite  differences.  The  choice 
of  a  higher-order  method  is  motivated  primarily  by  the  necessity  to  resolve  a  small-scale  phenomenon 
(backscattering)  at  the  background  of  the  forward  propagating  waves.  The  NLS  (2.10)  is  also  solved  by 
a  fourth-order  scheme;  it  is  natural  to  expect  that  this  will  leave  less  room  for  potential  purely  numerical 
discrepancies  between  the  two  techniques  and  as  such,  will  allow  for  a  more  accurate  comparison.  Besides, 
it  is  generally  known  that  higher-order  methods  provide  for  a  better  resolution  of  waves. 

3.1.  Numerical  Integration  of  the  NLH.  Our  numerical  method  for  solving  the  NLH  is  delineated  in 
[7];  here  we  only  outline  its  key  elements.  We  use  a  conventional  fourth-order  central-difference  discretization 
of  the  Laplacian;  in  so  doing  the  stencil  is  five-node  wide  in  both  r  and  z  directions.  As  the  equation  is 
nonlinear,  we  implement  a  nested  iteration  scheme.  On  the  outer  loop,  we  freeze  the  nonlinearity,  i.e., 
consider  the  coefficient  k 2  of  (2.1)  as  a  given  function  of  r  and  z7  which  is  actually  obtained  by  taking 
\E\ 4  from  the  previous  iteration.  This  way  we  arrive  at  a  linear  equation  with  variable  coefficients.  The 
latter  is  also  solved  by  iterations  on  the  inner  loop  of  the  nested  scheme.  Here,  we  leave  the  entire  varying 
part  of  the  equation,  which  is  proportional  to  e,  on  the  lower  level,  and  on  the  upper  level  need  to  invert 
only  the  constant-coefficient  linear  Helmholtz  operator  A  +  Formally,  our  iteration  scheme  resembles 
the  fixed-point  approach,  however,  no  rigorous  convergence  theory  is  available  yet,  and  the  convergence 
has  to  be  assessed  experimentally.  The  advantage  of  using  these  nested  iterations  is  that  first,  the  method 
eventually  reduces  to  the  repeated  solution  of  one  and  the  same  linear  constant  coefficient  equation  driven 
by  different  source  terms.  As  explained  below,  this  can  be  done  efficiently  on  the  discrete  level.  Second,  the 
radiation  boundary  conditions  and  the  two-way  ABCs  are  most  convenient  to  set  on  the  upper  time  level  of 
the  iteration  scheme  already  for  the  linear  constant-coefficient  operator. 

To  solve  the  linear-constant  coefficient  Helmholtz  equation  (discrete  counterpart  of  AE+k^E  =  g7  where 
g  is  the  right-hand  side  generated  on  the  previous  iteration)  we  first  separate  the  variables  by  implementing 
the  discrete  Fourier  transform  in  the  transverse  direction  r;  the  boundary  conditions  are  symmetry  at  r  =  0 
and  zero  Dirichlet  at  r  =  rmax.  This  yields  a  collection  of  fourth-order  one-dimensional  finite-difference 
equations  parameterized  by  the  dual  Fourier  variable;  each  of  the  latter  needs  to  be  solved  independently. 
The  two-way  and  radiation  ABCs  at  z  =  0  and  z  =  zIimX7  respectively,  are  set  in  the  Fourier  space  as  well, 
i.e.,  separately  for  each  of  the  aforementioned  one-dimensional  equations.  This  is  done  by  first  identifying 
the  linearly-independent  eigen-modes  for  the  homogeneous  version  of  each  one-dimensional  equation.  It 
is  important  to  note  that  even  though  the  original  differential  equation  is  of  the  second  order,  we  are 
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using  its  fourth-order  approximation  and  as  such,  each  homogeneous  discrete  one-dimensional  equation 
has  four  linearly  independent  solutions.  One  pair  of  the  latter  approximates  the  genuine  modes  of  the 
differential  equation,  those  may  be  either  traveling  or  evanescent  waves  depending  on  the  value  of  the 
dual  Fourier  variable.  The  other  pair  is  a  pure  numerical  artifact,  these  waves  are  always  evanescent,  but 
their  presence  implies  that  every  discrete  equation  requires  two  more  boundary  conditions  compared  to  the 
original  differential  equation.  The  radiation  boundary  conditions  are  constructed  by  requiring  that  on  the 
left  boundary  z  —  0  only  the  left-traveling  and/or  left-decaying  (evanescent)  waves  be  present  in  the  solution, 
and  on  the  right  boundary  z  =  £max  only  the  right-traveling  and/or  right-evanescent  waves  be  present  in  the 
solution.  The  selection  is  rendered  by  the  so-called  one-way  discrete  Helmholtz  equations,  which  are  the  linear 
homogeneous  relations  that  define  the  span  of  all  appropriate  modes  for  each  boundary.  The  two-way  ABC 
that  also  prescribes  the  incoming  signal  at  z  =  0  is  constructed  on  the  basis  of  the  corresponding  radiation 
boundary  condition  by  substituting  the  right-traveling  incoming  wave  into  the  one-way-to-the-left  Helmholtz 
equation  and  as  such  creating  the  inhomogeneity  of  a  particular  form,  see  [7].  Simple  considerations  based 
on  the  linear  superposition  principle  and  uniqueness  guarantee  that  the  resulting  nonhomogeneous  relation 
will  correctly  specify  the  incoming  signal  at  z  =  0  and  still  ensure  the  reflectionless  propagation  of  all  the 
outgoing  (i.e.,  left-traveling)  waves  through  z  =  0  toward  z  =  —  oo. 

As  concerns  the  computational  complexity  of  the  resulting  algorithm,  if  we  introduce  the  grid  dimensions 
Nr  and  NZ7  then  the  cost  of  both  the  direct  and  inverse  FFT  will  be  0(NzNr  In  Nr),  The  cost  of  solving 
each  of  the  Nr  one-dimensional  systems  will  be  linear  with  respect  to  Nz.  Indeed,  in  the  course  of  iterations 
each  of  these  systems  needs  to  be  solved  many  times  for  different  right-hand  sides.  Consequently,  the  sparse 
LU  decomposition  can  be  performed  only  once  ahead  of  time,  and  the  cost  of  each  backward  substitution  is 
linear.  Altogether,  the  complexity  of  each  iteration  is  still  0(NzNr  In  Nr). 

3,2,  Numerical  Integration  of  the  NLS.  The  NLS  (2.10)  is  discretized  in  the  r  direction  using 
standard  fourth-order  central  differences.  It  is  integrated  in  z  with  a  four-stage  Runge-Kutta  method  starting 
with  the  initial  data  Emc  (r,  0) .  The  boundary  condition  at  the  remote  lateral  boundary  r  —  rrnax  is  zero 
Dirichlet,  as  in  the  case  of  the  NLH. 

4.  Computational  Results.  In  accordance  with  the  discussion  of  Section  2  we  have  designed  a  set  of 
numerical  experiments  aimed  at  achieving  two  objectives:  (I)  Validate  the  computational  algorithm  of  [7]  for 
solving  the  NLH  through  a  comparison  with  numerical  simulations  of  the  NLS  model,  and  (II)  corroborate 
that  backscattering  effects  captured  by  the  NLH  model  scale  quadratically  with  the  nonparaxiality  parameter 
/,  as  suggested  by  the  analysis  of  Sections  2.3-2. 4.  Regarding  the  first  objective,  let  us  note  that  previously 
we  have  tested  the  numerical  algorithm  of  [7]  in  the  nonlinear  regime  using  grid  convergence,  but  never 
compared  it  with  any  other  algorithm  for  computing  the  propagation  of  waves  in  Kerr  media.  As  concerns 
the  second  objective,  it  amounts  to  the  accurate  numerical  computation  of  backscattering  in  nonlinear  self- 
focusing,  and  we  are  currently  unaware  of  any  previous  technique  with  similar  capabilities. 

To  be  able  to  conduct  an  accurate  comparison  of  the  numerical  predictions  obtained  with  the  NLH  and 
NLS  models  (equations  (2.9)  and  (2.10),  respectively),  we  have  chosen  a  regime  with  the  input  power  below 
critical,  for  which  the  solution  of  the  NLS  does  not  develop  singularities.  We  take  ko  =  8  and  e  =  0.04,  which 
corresponds  to  74%  of  the  critical  power  Nc,  see  [7, 11].  The  incoming  beam  profile  is  ^(r,  0)  =  e“(r/r°^, 
such  that  the  beam  width  r*o  is  much  less  than  rmax.  To  allow  for  the  variation  of  the  nonparaxiality 
parameter  /  =  l/&oro?  we  vary  the  beam  width  r q  while  keeping  the  wavenumber  kg  and  the  quantity 
that  controls  the  fractional  critical  power  unchanged. 

In  Figure  4.1  we  show  the  on- axis  amplitude  profiles  for  the  NLH  and  NLS  numerical  solutions; 
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Nonparaxiality  parameter  f=1/8 


Nonparaxiality  parameter  f=  1/1 6 


Fig.  4.1.  On-axis  amplitudes  of  the  solutions  to  NLH:  \E(Q7z)\,  and  NLS:  (er^fe^)  1^\i/}{Otz)\f  plotted  vs.  z[2Lof- 


Figure  4.1  (left)  corresponds  to  /  =  1/8,  and  Figure  4.1  (right)  corresponds  to  /  =  1/16.  We  plot  the 
values  of  the  computed  solution  on  the  axis  of  symmetry  r  —  0  because  this  is  the  most  interesting  location 
in  the  domain  where  the  genuinely  nonlinear  phenomena  take  place.  A  clear  manifestation  of  nonlinear 
self-focusing  is  the  “bump,”  or  peak,  on  the  solution  curves  in  Figure  4.1,  whose  value  is  higher  than  that  of 
the  incoming  wave  ^inc(0, 0). 

It  is  easy  to  see  that  for  both  /  =  1/8  and  /  =  1/16  the  NLH  and  NLS  curves  in  Figure  4.1  are  close  to 
one  another.  As  the  NLS  is  a  well-established  model  that  has  regular  solutions  for  subcritical  initial  powers, 
we  conclude  that  our  numerical  algorithm  for  solving  the  NLH  [7]  [that  starts  the  iteration  process  with  the 
initial  guess  E  =  0]  indeed  converges  to  the  correct  solution. 

We  also  notice  that  the  discrepancy  between  the  NLH  and  NLS  curves  on  Figure  4.1  (left)  is  larger 
than  that  on  Figure  4.1  (right).  This  behavior  is  expected  according  to  the  analysis  of  Section  2.3,  because 
the  discrepancy  between  the  two  curves  is  due  to  nonparaxiality  and  backscattering.  In  particular,  these 
simulations  suggest  that  the  NLS  is  indeed  the  asymptotic  limit  of  the  NLH  as  /  — y  0.  Perhaps  the  most 
apparent  manifestation  of  the  presence  of  backscattering  in  the  solution  of  the  NLH  is  that  the  computed 
value  of  the  total  electric  field  at  z  =  0  differs  from  that  of  the  incoming  wave,  as  one  can  clearly  see  in 
Figure  4.2  where  we  zoom  in  on  the  two  curves  obtained  for  f  —  1/8.  The  small  ripples  in  the  NLH  solution 
may  also  be  evidence  of  backscattering  (see  Section  2.3). 

Next,  we  quantify  the  backscattering  effect  by  computing  a  series  of  solution  pairs  (NLS  and  NLH)  for 
additional  values  of  /.  In  Figure  4.3(left)  we  show  by  asterisks  on  the  log- log  scale  the  quantity 


(4.1) 


(e*gr§): 1/4  E(0,Q)  -^(0,0) 


(e^)1/4£(0, 0)  —  A(0, 0) 


(cf.  formula  (2.8)),  where  E  is  the  computed  solution  of  the  NLH  and  ^  is  the  solution  to  the  NLS  that 
satisfies  the  initial  condition  ^(0, 0)  =  1,  for  /  =  1/16,  1/12,  1/10,  and  1/8.  The  solid  line  on  Figure  4.3(left) 
that  fits  closely  the  computed  data  is  0.75  *  /2'05.  This  essentially  corroborates  that  the  magnitude  of 
backscattering  indeed  scales  quadratically  with  the  nonparaxiality  parameter  /,  as  predicted  in  Section  2.3 
and  further  discussed  in  Section  2.4. 

Let  us  also  note  that  to  evaluate  the  quantity  (4.1)  we  do  not  really  need  to  solve  the  NLS,  because  the 
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Nonparaxiality  parameter  f=1/8 


Fig.  4.2.  Zoom-in  on  Figure  4.1  (left)  in  the  area  of  the  nonlinear  self-focusing  peak. 


initial  profile  at  z  =  0  is  given.  To  compare  the  actual  computed  solutions  of  the  NLH  and  NLS,  we  plot  on 
Figure  4.3 (right)  the  quantity 


(4.2) 


(e&oro) 1/4  |£(0,  z) |  -  |^(0,  z) | 


for  /  =  1/8  and  /  =  1/16  (same  values  as  on  Figure  4.1)  as  a  function  of  the  normalized  propagation 
distance  for  z  >  0.  Although  the  curves  on  Figure  4.3(right)  are  oscillatory,  we  still  see  that  that  the 
difference  between  the  solutions  of  the  NLH  and  NLS  decreases  for  all  z  with  the  decrease  of  /. 


Fig.  4.3.  Left:  The  quantity  given  by  formula  (4-1)  as  a  function  of  f  compared  on  the  log-log  scale  with  the  approximation 
0.75-/2-05  (solid  line).  Right:  The  difference  (4.2)  between  the  two  solutions  as  a  function  of  z/2Luf  for  f  —  1/8  and  f  —  1/16. 


5.  Conclusions.  We  have  compared  numerically  the  solutions  to  the  nonlinear  Schrodinger  equation 
and  nonlinear  Helmholtz  equation,  both  of  which  model  the  propagation  of  time-harmonic  electromagnetic 
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waves  in  Kerr  media.  The  NLH  was  solved  using  a  new  fourth-order  method  supplemented  by  the  two-way 
artificial  boundary  conditions  that  guarantee  the  proper  behavior  of  the  waves  as  they  enter  and  leave  the 
computational  domain.  As  the  NLS  is  considered  an  established  model,  the  agreement  of  the  NLH  and  NLS 
simulations  provides  a  good  justification  that  the  NLH  algorithm  indeed  converges  to  the  correct  physical 
solution.  On  the  other  hand,  the  NLH  is  a  more  comprehensive  model  that,  unlike  the  NLS,  takes  into 
account  the  phenomenon  of  nonlinear  backscattering.  As  such,  we  attribute  the  small  discrepancies  that  do 
exist  between  the  NLH  and  NLS  solution  to  nonparaxial  and  backscattering  effects.  By  analyzing  several 
computational  variants  that  correspond  to  different  values  of  the  nonparaxiality  parameter  /  we  have  been 
able  to  corroborate  that  the  magnitude  of  the  backscattered  wave  indeed  scales  quadratically  with  this 
parameter,  according  to  the  theoretical  predictions.  To  the  best  of  our  knowledge,  this  is  the  first  study  ever 
that  allows  for  an  accurate  quantitative  estimation  of  backscattering  in  nonlinear  self-focusing. 
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